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Abstract- We present a framework for simulating signal propagation in geometric networks (i.e. 
networks that can be mapped to geometric graphs in some space) and for developing algorithms 
that estimate (i.e. map) the state and functional topology of complex dynamic geometric net- 
works. Within the framework we define the key features typically present in such networks and of 
particular relevance to biological cellular neural networks: Dynamics, signaling, observation, and 
control. The framework is particularly well-suited for estimating functional connectivity in celhilar 
neural networks from experimentally observable data, and has been implemented using graphics 
processing unit (GPU) high performance computing. Computationally, the framework can simulate 
cellular network signaling close to or faster than real time. We further propose a standard test set 
of networks to measure performance and compare different mapping algorithms. 
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1 Introduction 

Complex dynamic networks permeate many real world engineering and biological applications. The 
development of mathematical and computational tools for understanding and predicting network 
dynamics will be key to manipulating and interacting with such real world networks. Network 
theory is a subset of graph theory where the connections between vertices have a number value 
describing some attribute of that connection, such as for example bandwidth, flow rates, or a cost 
function. Complex networks are defined to have a non-standard topology, i.e. the functional links 
between nodes in the network, implying some structure in the connectivity pattern of the network 
beyond a simple lattice or complete random connectivity. Biological cellular neural networks are 
both complex and dynamic, meaning that the connection attribute between any two vertices may 
change with respect to time and, more importantly, individual vertices exhibit their own nonlinear 
signaling dynamics. Complex functional interactions of networks made up of large numbers of neu- 
rons and glia produce emergent systems-level phenomena such as consciousness and self-awareness, 
and are responsible for how neural information is represented and processed. Changes in the struc- 
ture of such networks presumably underlie the development of multidimensional central nervous 
system disorders. For example, hypersynchronous neuronal and glial activity in networks of neu- 
rons are associated with the paroxysmal depolarization shifts that underlie epilepsy [3H UIJ HSJ [12] . 
Ultimately, the physiologic behavior of a neural cell network is dependent on both its functional 
topology and the dynamics of individual cells. 

Within a complex dynamic network there are two topologies. A static, structural topology 
that describes all the possible connections within the network, and a dynamic, functional topology 
that establishes how a signal propagates through the static topology. Functional topologies are 
subsets of the structural topology and vary depending on the functional connectivity, internal 
dynamics of individual vertices, and the specific stimulus to the network. While this is the case for 
biological neural networks, where cells that are physically connected need not necessarily signal each 
other, in cellular neural circuits and networks structure and function influence each other and the 
states of cells and the connections between them may change with time as a function of plasticity 
mechanisms. However, structural changes in the physical connectivity of a cellular neural network 
leading to changes in the connectivity topology occur on a very different time scale then functional 
changes that can be influenced relatively quickly by plasticity mechanisms that produce changes in 
signaling efficacy between cells (i.e. changes in connectivity weights). While the observation of the 
structural network topology of cellular neural networks may be experimentally very challenging (and 
indeed is the focus of much intense research) , it is a relatively straightforward task. The observation 
of functional topologies in biological neural networks however poses additional experimental and 
theoretical challenges that need to be considered. Signaling events and resultant networks may 
be unique and be observable only once as a signal propagates through a network. The functional 
topology is dynamic and may change during observation. Noise and unknown external factors limit 
observability and reduce repeatability. These factors make the estimation of functional connectivity 
from observed activity a difficult task, though a critical one for systems neuroscience if we are to 
understand how dynamic functional signaling in the brain at the level of networks and circuits 
produces responses and behaviors in the organism. 

Current approaches for studying cellular neural networks can be roughly classified into three 
categories. The first and most popular amongst experimentalists are statistical methods that cor- 
relate the activities of two or more neurons in a network. This provides purely descriptive statistics 
about the behavior of cells. For the most part, statistical approaches make no underlying assunip- 
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tions about the cellular and systems dynamics that give rise to observed signals in a network of 
cells. Another way to study networks is through simulation of networks with known connectiv- 
ities and dynamic parameters in order to simulate real-world observed system level phenomena 
such as vision and audition. Using well established environments like NEURON or Genesis, many 
real-world phenomena have been described through simulation. However, dynamic parameters and 
functional connections are manually specified in simulation environments such as these in order to 
achieve results that mimic biological function, requiring the estimation of experimentally unobserv- 
able variables. The third category is in some ways the reverse process to simulation, where temporal 
data is used with appropriate models in order to estimate parameters. Within this third category, 
we introduce a modeling framework for using real-world data to map the functional topology of 
complex dynamic networks. While not a mapping algorithm or simulation environment, the frame- 
work formally defines key features of cellular neural network signaling and experimental constraints 
associated with observation and stimulus control, and can accommodate any appropriate model of 
intracellular dynamics. Alongside the definition of the framework, a test set of synthetic networks 
with known connectivities is provided to help the development of mapping algorithms by providing 
a common benchmark any such algorithm should be able to map. In a subsequent paper to this one 
we will introduce an approach that will estimate and map the functional topology of complex net- 
works with unknown connectivities given limited and often noisy observations that takes advantage 
of the results introduced here. 

The proposed framework has a number of unique properties that makes it particularly appli- 
cable to the constraints and experimental limitations imposed by real biological cellular neural 
networks. First, dynamic activity and signaling is modeled at the individual node (i.e. cell) scale. 
The dynamics of individual cells are modeled as state sets, with transition functions describing 
their evolution across discrete time steps. Cellular resolution was chosen because it represents 
the best compromise between observability, dynamics, and complexity. Large numbers of individ- 
ual cells can now be observed in parallel in functional neural networks using optical microscopy 
[M[ [SI HU [53] . Single cell neuronal dynamics are well understood and many models exist (see for 
example [T^ 143) ). while similar models of single cell astrocyte dynamics are beginning to emerge 
[ini [32l [35j [5l [28] . Attempting to go to a finer, sub-cellular compartmental resolution dramatically 
increases the complexity of the model, computational demand, and is generally not experimentally 
observable at a network level. Secondly, cells are located in physical space and their positions 
are easily determinable during experimental observation. When connected cellular networks form 
geometric networks. Thirdly, the effect of a signal on a target cell is defined as a state change in 
the target cell in response to the influence of a source cell that connects to it. That influence is not 
instantaneous, and is delayed by the physical distance between cells and the speed of transmission. 
Signals are modulated in strength by functional weights, which establish the magnitude of the in- 
fluence. Fourthly, to more realistically simulate experimental conditions and measurements, noise 
can be added to multiple levels within the framework, from parameters to state and observation 
variables. Finally, experimental user-defined controls at the individual cell level are defined within 
the framework. Controls should be designed to make observations more informative of the network 
dynamics, but should not change the underlying parameters and connectivities. The framework is 
described in detail in section[2] The results section (section[3| shows how single cell dynamic models 
are integrated within the framework (3.1), and how network connectivity is established from indi- 
vidual cells (3.2 1. We also describe how the framework accommodates plasticity mechanisms (3.2.3 1 
and experimental observability associated with optical calcium imaging (3.3). Section 3.4 discusses 
the practical implementation of the framework using high performance graphical processing unit 
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(GPU) computing. 

In section 3.5 we use the framework to propose a standard set of benchmark test networks 
of varying sizes and topologies to evaluate and compare different network mapping algorithms. 
Mapping algorithms would have access to simulated observable data (i.e. simulated experimental 
data) generated by the framework as a function of chosen test networks and be required to derive 
the unobservable parameters and functional network connectivity. The concept of a standardized 
test to gauge the effectiveness of an algorithm is not new, especially for optimization algorithms. 
For example, in the field of nonlinear programming and optimization a standard benchmark set was 
established in a landmark collection of test problems |23j that are used for testing any nonlinear 
optimization algorithm. Test collections have grown and developed into problem environments, 
providing the underlying problem code to be used directly by the optimizers [HI HI])- By providing 
a set of problems with known solutions, algorithm developers have a standard by which to measure 
solution accuracies, convergence rates, computation times, and suitability to different problem 
types. We propose that a similar test set for algorithms designed to identify and map functional 
cellular neural networks and circuits will be just as useful. To address this, we have developed 
computer code that generates observable data from a known network and connectivity. The code 
encompasses all the elements of the framework, runs in real time for all the test networks, and 
is designed for parallel computation, and can therefore be used as a starting point for mapping 
algorithms. 



2 A Framework for Dynamics, Signaling, Control and Ob- 
servation in Geometric Networks 

We develop the proposed framework using standard graph theoretic and set theoretic concepts and 
terminology. In the most general sense, a network is a type of graph. A graph is defined as an 
ordered pair of finite disjoint sets (■;;, E) such that v is the set of J vertices of G and E is the set of 
edges of G, i.e. v — v{G) is the vertex set of G while E = E{G) is the edge set of G. An edge Cij 
is defined if there is a directed connection from vertex i to vertex j. Geometric graphs are graphs 
where the relative positions of vertices are assigned coordinates in some geometric space. While 
this is the most generic description of a graph, dynamic geometric networks as we use the term here 
are more specialized cases of generalized geometric graphs defined as follows. Vertices in a network 
have two attributes, a known and static position in physical Cartesian space denoted by Xj for a 
given vertex j and a time- variant state set yj{t) of Kj state variables: 

yj{t) = {yi,At),y2,jit), . . .,yK,,j{t)} (la) 

such that formally 

yj{t) = {yk,j{t) : k G N,k < Kj} for any given vertex j (lb) 

Next, for all vertices i other than j, let the set Yj{t) be the union of all i, i.e. the collection of 
states of all vertices in the network excluding vertex j, weighted and delayed relative to vertex j, in 
the sense that every vertex i has the potential to pass information (e.g. a signal) to vertex j with 
varying amounts of 'influence' as determined by a collection of weights that modulate any directed 
edges from i to j. Furthermore, such information will be delayed by some finite time as a function 
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of the geometric position of vertex i in the network relative to j and the finite speed of information 
propagation. We define 

Yj{t) = UigN;i<J;j^jf^ij(t) • yi{t " Tij)} (2) 

where with out loss of generality we define 

it)] (3) 

and restrict [Tb| for vertex i with temporal delays as vector sets, i.e. 

yi{t) = [UkAt) ■■ k eN,k < Ki] for any given vertex i ^ j (4) 

The delays are non-negative values representing the delay of information passing from i to j. In 
all cases, here and below we adopt the convention that indexing subscripts given by 'ij' enumerate 
the variable that uses the subscript as linking vertex pairs i and j. 

We then define a transition function Hj(-) with parameter set Qj that describes the temporal 
progression or evolution of yj{t) in discrete time increments At: 

y,(< + Ai) =H,(y,(t),Y,(i),u,(t),e,) (5) 

where Hj(-) is given by 

Hj = UkeN;k<KjHkj{yj{t),Yj{t),Uk,j{t),Qk,j) (6) 
Uj(t) is a user control or experimental input. 

= UfeeN;fc</fjMfcj(i) (7) 

and Qj is parameter set 

Okj = {Oik :l,keN;l< Lk] k < K,) (8a) 

9j = UfegN;fe<A:j6fej = ^k=i^k,j for any given vertex j (8b) 

and 0,7 = Uj6N;j<jej = u/^iOj (8c) 

Lk is the number of parameters for a given state variable, Kj is the number of state variables for 
a given vertex j, and J represents the size of the network (i.e. the total number of vertices). Note 
that the functions comprising the set Hj(-), each advance their respective variables in time: 

yij(i + At) = i7ij(?;ij(t),Yj(t),ui,,,eij) 

y2.j(t + At) = i72,,(2;2j(i),Y,(t),U2,,,e2j) 

J- (t + At) = Hk^ ,i (2/i^^ ,j (t) , (t ) , UK^ J ,Qk,,j) 

Similarly, we define a function Gij(-) that describes the time course of the weighing sets ^ij{t) 
with parameter sets Aij as follows: 

%(t + At) = Gy(%(t),yj(t),y,(t-T,j),A,,) (9) 

where analogous with|6]Gij(-) is given by 

Gij = UkeN;k<KGk.j{^ij{t),yj{t),yi{t — Tij),-^-k,ij) 
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with parameters 

^kaJ = {A;,fe ■.l,ken-l<Lk]k< K} (10a) 
Aij = i^ken-M<KJ^k,ij for any given vertex pair ij (10b) 
The delays between vertex pairs r^- are defined as functions of the positions of the two vertices: 

r,, =i^(x„x„r,,); D{-)>Q (11) 

where the set Tij is the set of parameters of the non- negative function _D, specific to the pair ij. 

Formally, the temporal evolutions of {t) and ^lij (t) are continuous and expressed as a discrete 
delay differential equations with delays Tij for all vertices connecting to vertex j, so that the 
continuous forms of equations [5] and [9] are 

^=H,(y,(t),Y,(i),u,(t),e,) (12) 

and 

G,, (f],, (t), y, (t), y,(t - r,,), A,,) (13) 

In the limit as At — > [s] and [o] can be written as yj{t + dt) and flij{t + dt) and equations 12 anc 13 
apply. However, from a practical experimental perspective time measurements will always be finite 
and the discrete forms need be considered. As such, in this paper we do not pursue further the 



interesting theoretical implications of the continuous forms given by 12 anc |13| 

Finally, we define an observation set 2j{t) composed of M variables, that operates directly on 
the state set yj{t): 

z.W=F,(y,(0,$,) (14) 
The observation function Fj is vector valued: 

Fj{t) = UmeN-m<MF,n,j (yj(t), $mj) (15) 

with parameter set $j given by 

= {<l}i,m -.1,1 en-J < L^;m< M} (16a) 

*j = Ur,ieN-m<M^m,] for any given vertex j (16b) 

The framework presented here is general, as it allows for communication between any two 
state variables between any two vertices. Transition functions and their parameters are defined 
specific to vertex Hj(-) or communication between vertex pairs Gy (•). This produces a large set of 
functions and parameters, though in practice one or two different functions are applied to all cells 
or combinations. The weighing set can operate on all state variables of the connecting vertex i 
into target j, though usually one state of i is transmitted to one state in j. In the next section, we 
will describe how several dynamics and communication models used in cellular networks fit within 
this framework to reproduce observable quantities similar to experimentally measured data. 
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3 Results 

The framework can accommodate essentially all models of both neuronal and astrocytic dynamics. 
Independent of the specifics of any single cell model chosen, the framework provides a compact 
mathematical structure that quantitatively describes signaling and information propagation and 
flow in geometrically defined networks. The geometry and physical connectivity topology of the 
network can be simulated (e.g. random, scale free, or small world) or measured from experimen- 
tal data such using methods such as optical imaging. Regardless of how one choses to set up the 
network, the framework provides a description of information flow through the network given knowl- 
edge of temporal signaling delays and chosen single cell models, or can be used to identify and map 
unknown functional connectivities and parameters in real neural circuits and networks. In all cases, 
the framework is able to provide an estimate of the complete description of the functional network 
and the interaction between all observable and hidden state variables and parameters. Figure [T] 
illustrates a simple five vertex example that summarizes everything that is needed to describe the 
functional dynamics of information flow through the network. Figure [2] provides a specific example 
of the network from figure [T] using a Hodgkin- Huxley model and simulating one second worth of 
data. Note how the framework provides experimentally measurable variables (calcium and mem- 
brane voltage) for every cell in the network in the temporal sequence dictated by the geometry and 
connectivity of the network. 




Figure 1: A five- vertex dynamic network. Each vertex j has a position in physical Cartesian space 
denoted by vector Xj, and a dynamic state set Yjit). A vertex's dynamic state varies in discrete 
time steps, and is influenced by its own previous state and the states of other vertices connecting 
into it, with a delay r and a functional connection weight n. The time delays between vertices are 
a function of their positions in space. The magnitude of the connection weights, ^lij, are estimated 
based from vertices' known positions and the observed dynamics. 



3.1 Individual Cell Dynamics 

In this section we discuss how neuronal models of single cell dynamics, synaptic connections, plas- 
ticity, and observation flt within the framework. We begin by showing how to construct the state 
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Figure 2: Signaling dynamics of the network from figure [T] A Hodgkin-Huxley model is used for 
single-cell dynamics in a one second simulation. The delays (written in blue) are in milliseconds 
and are based on the cartesian distances between cells. Functional weights, shown as relative 
strengths in green, are chosen arbitrarily for the purposes of this example. The right panel shows 
the network dynamics for each of the five cells following a sequential pulse stimulation (red bar) at 
each cell. Experimentally observable calcium traces are shown as solid black lines on an arbitrary 
vertical scale. Voltage, a hidden variable, is shown as a gray line and constitute the neurons' action 
potentials. 



transition of an isolated (unconnected) vertex, and then build the full transition function by con- 
necting multiple vertices into a network. In an isolated vertex case with no incoming connections, 
the state transition reduces to 

y(i + A<) = H(y(t),u(t),e) (17) 

This generic form encompasses neuronal models described in differential equation form, as well as 
those with a state reset based on some threshold value. Most neuronal models are expressed in 
differential form as 

^^h{y{t),n{t),e) (18) 

Converting [18] into the state transition form given by [17] is a matter of numerical integration with 
an integration method of choice. Using Euler's method, for example, the state transition function 
Hj of the system in 18 becomes 

H,{y,{t),u,{t),e,)=y,{e,t) + At-h{y,{t),u,{t),ej) (19) 

Here we used the Euler method of integration for its simplicity and clarity, but other, but more 
complex integration methods like trapezoidal or Runge-Kutta can also be used to generate the next 
time step from the current step. 

As an example, consider the simple Fitzhugh-Nagumo oscillator used to model neurons. In its 
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differential form the model is given by the pair of equations 

aV{t) - hV{tf - cW{t) + S{t) + U{t) (20) 
e{V{t) + f- gWit)) 



dt 

dW{t) 



dt 

The state set y{t) is comprised of two state variables yj{t) ~ {Vj (t) , Wj (t)} . The experimental 
control set is composed of only one variable, affecting the V{t) state variable, so Uj{t) = {Uj{t)}. 
This system can be expressed in state transition form as Hj = {Hyj, Hwj}, with 

Hv,j (t) = Vj (t) + At ■ {aVj (<) - bVj {tf - cWj (t) + Uj (t)) (21) 
Hw^j{t) = Wj{t) + At ■ {e{Vj{t) + f- gW,{t))) 

Note that we index the state transition functions based on the state variables they operate on; for 
example, Hyit) advances V{t). The parameter set for this system is composed of the parameters 
for each of the state transition equations in [2l] 

©Vj = {a,b,c} 
Qw.j = {ej,g} 

e, = ei., U 62,, (22) 

This system has two state variables and six parameters. 

Another class of neuronal models are those with a hard reset. These models are also described 
in differential equation form, but contain a hard reset when a state variable reaches a certain value. 
As an example, consider the Izhikevitch simple model, written in differential form as 



dV{t) _ 1 

dt ~ C 
dW{t) 



k(V{t) ~ Vr){V{t) - Vt) - W{t) + S{t) + C/(t)J (23) 
a{h{V{t) - Vr) - W{t)) 

if V[t) > Vpeak (spike event) 



dt 

V{t+) ^ c 
W{t+) W{t) 

This model consists of a voltage and amplifying currents {V{t)) and a resonant gating variable 
{W{t)). The system has up to 9 parameters, and resets both state variables when a certain voltage 
threshold {Vpeak parameter) is surpassed. Mapped onto our framework, the model and its transition 
functions are 



H^^^ ^ ) V, it) + At • i [k{V, {t) - Vr){V, (t) ~Vt)~ W, {t) + U, {t)\ if V, {t) < Vpeak ^^4) 

Vreset Otherwise 

^ W, (t) + At ■ [r{b{Vj it) - Vr) - W, (t))] if Vj < Vpeak 

^'^ \ Wj{t)+d otherwise 

Here, the parameter set is Qj — {C, k,Vr,Vt, Vreset, Vpeak} U {r, 6, c?, V^eafe}, and just like the 
Fitzhugh-Nagumo model, yj{t) = {Vj{t),Wj{t)}, Uj{t) = {Uj{t)} and = {Hv,j,Hw.j}- By 
changing the parameter values of the individual models, different classes of neurons can be simulated 
with the same transition function. 
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Similarly, any model can be accommodated and fit into the framework, from the simplest to 
the most complex. Traditionally, all neuronal models have membrane voltage as a state variable 
and propagate a discrete signal in the form of an action potential when the membrane voltage 
rises past some threshold value at the axon hillock in response to depolarizing and hyperpolarizing 
currents in dendrites mediated by spatial and temporal summation of presynaptic currents. The 
simplest neuronal model, the leaky integrate and fire (LIF) has voltage as a single state variable 
that decays to a target value and is perturbed by incoming currents. If the voltage rises past a 
threshold value it is reset at the next time step to a reset value. One of the most complex and 
realistic single-cell models is the Hodgkin-Huxley (HH) model which relies on four state variables 
{v, m^rijh} to describe the dynamics responsible for the generation of action potentials. The number 
of parameters increases with the number of state variables, from 4 in the LIF model to 22 for the 
HH model. Additionally, the required time step is shorter for HH models, being on the order of 0.03 
milliseconds compared to roughly 5 milliseconds for the LIF model. The increased number of state 
variables and parameters along with shorter time steps puts a significant computational burden on 
any simulation or mapping algorithm. The question of which model and how much complexity is 
required to best describe real-world data is not trivial and depends on the purpose and intent of 
the modeling. 

Although considerably more limited than the number of existing and studied neuronal dynamic 
models, there are a few astrocyte dynamic models that emphasize differing aspects and processes of 
astrocyte signaling [JDl 13211351 13 UHl ISO]- Astrocyte models are expressed in differential forms similar 
to equation[T8j Further research into astrocytic models is important though because astrocytes have 
been shown to play a direct role in the bidirectional communication between themselves and neurons 
via intracellular calcium transients and intercellular calcium waves under controlled experimental 
conditions [20l [1] HH [121 EHl [37] and more recently physiologic conditions in the neural retina [27] 
and cerebellum [25]. Pathophysiologically intercellular calcium waves in astrocytes independent of 
neuronal hyperactivity have recently been shown to occur spontaneously in vivo in the APP/PSl 
transgenic mouse model of Alzheimers disease [35]; and amyloid beta has been shown to be sufficient 
to trigger complex temporally delayed intercellular calcium waves in isolated astrocyte networks 

The state transition framework handles all single-cell dynamic models, as well as heterogenous 
systems of different cell types, either by different parameter sets or state transition function sets 
or both. Whatever the dynamics of individual neurons or astrocytes, all perform the same general 
task whereby processes and inputs generate outputs to other cells in a connected network. 

3.2 Cellular Network Signaling 

There are three components to cellular signaling: how long it takes for information from one cell to 
reach another, what are the effects of one cell on another and how do those effects change through 
time given the relative dynamics of the two connected cells. In this section we describe how signaling 
delays, functional connectivity and plasticity are accomodated by the framework. 

3.2.1 Signaling Delays 

Neurons and astrocytes form signaling networks that pass and process information between cells, 
and the state transition function must be extended to include signal propagation between func- 
tionally connected cells. In biological cellular networks, signal propagation occurs at a finite and 
relatively slow speed (i.e. compared to electronic circuit networks). Information in cell networks 
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propagates on the order of microns per second for astrocytes to meters per second for myelinated 
axons in neurons. Thus the influence of the dynamics of one cell is felt by another cell after some 



delay t^ . While the general form for the delay is given by 11 above, the simplest form it can take 
is when signaling is a geodesic between vertices (or between the centers of cell bodies or centers 
of other cellular compartments as needed in morphologic models) and the transmission speed s is 
constant: 

= (25) 



Tij = D{T,:xii,Tij), = {s} (c./jTT] above) 

Here the delay is simply the Euclidean distance between the cell centers divided by the trans- 
mission speed. For a diffusive network, as is the case with astrocytes, delays are proportional to 
the square of the distance between vertices. A more complex delay function may take into account 
knowledge about the particular physiology of the network, curved paths between cells, non-uniform 
speeds, etc. The dependency of the framework on the delays is critical to its ability to describe how 
and when information within the network is processed, ultimately to a significant degree dictating 
the intercellular dynamics of the overall neural circuit or network. 

Figure [3] illustrates the dependency of network dynamics on signaling speed and delay times in a 
100-vertex three dimensional network by varying the intercellular signaling speed, with everything 
else, including its geometry (i.e. its physical connectivity), the functional connectivity and input 
stimulus, remaining the same. We stimulated all cells with 500 ms of depolarizing current. The 
delays are inversely proportional to the signal propagation speed. We illustrate the effects of three 
signaling speeds, 2, 20, and 200 pixels/ms. At the lowest signaling speed, 2 pixels/ms, a low- 
frequency periodic activity was produced that qualitatively resembles a central pattern generator. 
At a speed of 20 pixels/ms fewer cells exhibited low frequency oscillations, and signaling became 
more sporadic. For both the 2 and 20 pixel/ms propagation speeds however, signaling continued 
past the period of stimulation. At a propagation speed of 200 pixels/ms however, there was no 
signaling past the stimulus period. In fully recurrent networks such as the one illustrated here, 
delays serve as a form of signal storage, essentially giving cells time to recover from a refractory 
period between activations, which in turn maintains recurrent signaling propagation well beyond 
an initial stimulus. For some appropriate range of signaling speeds, and therefore delay times, this 
recurrent signaling can settle into a repeatable pattern. If however, the signaling speeds are too 
fast, incoming signaling from upstream cells never have an opportunity to activate downstream 
cells because they are still refractory and do not respond. This leads to signaling in the network 
quickly dying away and not being sustained without it being driven by an external stimulus, as 
is the case with speeds of 200 pixels/ms in this example. A full discussion of the dependency of 
the network dynamics on the variables that govern it is very involved and beyond the scope of 
this paper. However, this example serves to illustrate that along with a network's connectivity 
topology and individual vertex dynamics, signaling delays play a crucial role in its overall response 
and dynamics and must be part of any network simulation or modeling framework that attempts 
to capture the inherent behaviors of neurobiological networks. 

3.2.2 Fiinctional Connectivity 

The other component of signaling is the functional connectivity of the network, or how the state 
of one vertex influences the state of another. In equation [2] the set collects all the states of all 



11 



Buibas and Silva 



Simulating and estimating dynamic geometric networks 



2 pixels/ms 



500 pixels 



B 




20 pixels/ms 



20 

"5 40 



o 60 



80 
100 



200 pixels/ms 



20 



40 60 80 100 
Target cell 



100 
80 

60 f 

40 i 

o 

20 

100 
80 
60 
40 

20 

100 

80 

60 

40 

20 



500 1000 1500 2000 2500 3000 3500 4000 
Time (ms) 



Figure 3: Effects of signaling speed on network dynamics. The network with spatial locations and 
physical connections shown in panel A, is assigned random weights uniformly distributed between 
-1 and 1 on each physical edge, panel B. An Izhikevitch simple model of bursting neurons was used 
to model the individual vertex dynamics. C. By varying the speed of signal propagation, the delay 
distributions are scaled, having a substantial impact on the spike dynamics (see text). 



12 



Buibas and Silva 



Simulating and estimating dynamic geometric networks 



the vertices in the network except j, delayed by a time value relative to j. When Yj is passed 
into the transition function Hj, the information contained in every state of every other vertex is 
made available to all state variables in vertex j. Within the framework, this is the broadest possible 
scope of connectivity, though in practice typically information from one variable affects one or more 
variables in another vertex. 

As an example, consider a neuron and its pre synaptic (chemical) inputs that induce post 
synaptic currents. The post synaptic current (PSC) is modeled with additional state variables in 
the neuron's state vector yj{t). There are two PSC models in wide use which we have tested within 
the framework, although it is in no way limited to these two examples. In general, a pre synaptic 
neuron causes a PSC in a post-synaptic neuron that ultimately affects membrane voltage. Modeling 
the effect of an incoming signal on a target cell is key to establishing connectivity based on the 
observed cell dynamics. Post synaptic current is represented as a function of the form 

sit) = g,nax ■ r{t) ■ {V{t) - Erev) 

where the resultant signals sj (i) are summed and passed to the voltage state variable in connected 
vertex i. gmax is the maximum conductance, and in the case of a specific synapse, can be expressed 
as a product of the maximum allowable conductance and the (instantaneous) functional connection 
weight Uij. r(t) describes the time course of the current, and is generally one of two forms, depending 
on the neuron type and neurotransmitter release |151 116| : 

r{t) ~ e^°* simple exponential 

r(t) = at ■ e""* a-function 

where a is the time decay constant. In both cases, time starts at the moment of activation, in this 
case the time of the activation of the pre synaptic cell plus the delay to the post synaptic cell, Tij. 
For either case, the expressions for r{t) can be written as linear differential systems, with the spike 
as the impulse. For the simple exponential the differential equation is 

drAt) , , 

rj{t+) = rj{t) + LOij upon arrival of spike event from vertex i (26) 
For the a-function, another state variable p{t) is used: 

= -ah,{l) - 2ap,(t) 

Pj(t+) — Pj{t) + ujij upon arrival of spike event from vertex i (27) 

Shifting the arrival delays is simply a matter of shifting the spike detection function of the pre 
synaptic neuron, so the arrival time from vertex i to vertex j is effectively the time shifted function 
of the voltage of i : Vi{t — Tij). The synaptic current is a decreasing exponential with rate a, 
incremented by a weight value ujij upon arrival of a spike occuring Tij time units ago at another 
cell i. 

For example, extending the state set for the Izhikevitch model given by[24]for an arbitrary vertex 
to include an a-function PSC model of post synaptic currents produces yj{t) — {V, (i), Wj{t), rj{t),pj{t)} 
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and Hj — {Hvj, Hw,j, Hrj, Hpj} where, 

Hy =1 ^1 (^) +^^-h [^C^J- (*) - ^r){V, (t) -Vt)- Wj (t) + r{t) + Uj (t)] if V, {t) < Vp,ak 
\ Vreset Otherwise 

(28) 

^( W, (t) + At ■ [r{b{V, (t) - Vr) - W, (t))] if Vj < Vp,ak 
\ Wj{t) + d otherwise 

Hrj ^ rj{t) + At ■ pj{t) (29) 
Hp,j = Pjit) + At- I -a^Tjit) - 2apj{t) + ^ LOkjSvA{Vk{t - Tfe,)) ) (30) 



The spike detection function is given by 



'^'nvim^W IZ^'- (31) 

Here, the spike inputs from other neurons are passed into the state variable p{t) with transition 
function Hpj, integrated by variable r{t), and passed into the voltage variable V{t). The functional 
weights only operate on the voltage spikes from other neurons, so only information from the voltage 
state variable V{t) is passed to other neurons. A similar set can be constructed with the PSC 
models in [26l ■ 

3.2.3 Plasticity 

Up to this point we have described models with fixed functional connectivity strengths. This 
assumption is valid for networks observed over short periods where connection strengths can be 
assumed to be constant for the purposes of mapping or simulating since the plasticity mechanisms 
that modulate connection strengths operate on longer time scales. If the weights change as a 
function of the activities of the cells it connects, the framework can be used to modulate connective 
strengths (equations|9]and 10 ). Just as with single cell dynamics and network connectivity, there are 



many models of plasticity and we will not attempt to list or review them all here. Rather, we describe 
how a simple spike-time dependent plasticity model reviewed by Bi and Poo in [7] that incorporate 
long-term potentiation (LTP) and long-term depression (LTD) can be easily implemented within 
the proposed framework. 

In equation [9] we described the functional strength or weight Vlij with transition function Gy- 
analagous to the state set for individual neurons. When the weight is constant, the set contains 
only one variable so ^lij — {wjj} and there is no temporal change in ujij and no transition function 
or parameters. However, if the weight is modulated by the activitiy of cell i on j (the connection is 
directional, so uiij modulates information flowing from i to j), then the states yi(i — r^j) and yj(i) 
affect how VLij (t) changes in time. 

The neuronal plasticity model in i7j describes strengthening and weakening of synaptic conduc- 
tance based on the timing of spikes in pre and post synaptic neuronal spiking. If a post synaptic 
spike occurs immediately after a pre synaptic spike, the synaptic conductance is increased and the 
functional connection is effectively strengthened. If the post synaptic neuron spikes before the pre 
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synaptic neurons, the connection is weakened. Other conditions hke spike coincidence or long times 
between the spikes of pre and post synaptic neurons have no affect on the synaptic conductance. 

To incorporate this model into our framework we first augment the set of cell states with another 
state variable Sj{t), that stores the time from the last spike. This state variable has transition 
function Hsj: 

H,^j - (1 - spd(y,(t)) • {sj{t) + At) (32) 

The model in [7] uses two exponential curves to describe changes in synaptic conductance based 
on the spike times. By defining the function qij{t) as 

Qijit) = {si{t - Tij) - Sj{t)) ■ max(spd(Vi(i - nj) , spd{Vj (t)) 

we can reconstruct the plasticity model within the framework as 

= a ■ sign(<z,,(t)) . exp (-/3|g.,(t)|) 

Thus, the transition function for dynamic synaptic weights is 

G,y = + At-{a- sign(g,j(t)) • exp (33) 

The parameter set for connection ij is A.^- = {a, /?}. This way, a static weight ujij can be converted 
into a dynamical one ujij{t), with behaviors governed by any arbitrary plasticity model and its 
associated parameter set A^j. It is important to note that while here we describe only one scalar 
weight between two vertices, the framework as defined can accomodate as many weights as there 
are state variables for a given vertex, thereby describing different classes of intercellular signaling 
between vertex pairs. For example, in networks of neurons signaling may occur via gap junctional 
mediated electrical synapses in addition to chemical synapses, while in astrocyte networks inter- 
cellular signaling is typically mediated by diffusional processes (e.g. vesicularly released adenosine 
triphosphate, ATP). 



3.3 Experimental Observability Through Calcium Observation 

Typically only one or a few of the state variables in the state vector are available for observation 
(i.e. are experimentally measurable). This is certainly the case with cellular networks, especially in 
neurons where voltage is measured as an indicator of signaling activity. But simultaneous voltage 
measurements are difficult for networks of many of neurons where the geometry of the network may 
be important to the analysis or interpretation of the data. While high density planar multi-electrode 
arrays can record from a few hundred cells at once, it is typically not possible to correlate recorded 
activity with the native geometry of the network. (There is one notable exception to this that we are 
aware of, which is the ganglion cell monolayer in the peripheral retina. Chichilnisky and colleagues 
are able to computationally infer the geometry and functional activity of these retinal ganglion cells 
due to their unique planar arrangement- see [34] EHl [10]. However, in the brain and even in the 
retina where mutlilayered ganglion cells recieve incoming macular inputs it is not possible to do with 
electrode arrays.) Given the challenges associated with direct electrophysiological measurements of 
large neuronal ensembles, calcium fluorescence imaging has been used as an indirect measure of the 
dynamics of large neuronal netwokrs. The dynamics of calcium, often influenced by voltage spiking, 
are modeled as state variables and associated transition functions added to the state set. There 
are a few models describing the time course of calcium as it is driven by changes in membrane 
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voltage, and attempts to develop more refined and robust calcium dynamics models specifically for 
the study neural microcircuits and networks is a very active area of research. Our intent here is 
not to describe on-going efforts or the state of the art but simply to illustrate the integration of 
one such model within our framework. 

The simplest model of calcium dynamics can be expressed as a linear system, with a spike input 

T,^=-lc(t)+spd(y(0) (34) 
at Tc 

This model is integrated into the framework by the addition of another state variable c(i) describing 
the cytosolic calcium concentration to the state set, with transition function: 



Hc^j = c,(t) + At (^-^c,(t) + spd{Vj{t))^ 



This model has one parameter Tc which describes the removal rate of calcium, after an input caused 
by a spike. This is the simplest model used for calcium dynamics based on neuronal spiking and 
is often used to extract spikes from calcium |45| . More complex and non-linear models of calcium 
dynamics have been developed and these can also readily be integrated within the framework . 

Finally, there is the issue of observable variables. In equation [14] we defined an observable set 
Zj{t) as some function of the current state set yj{t). When using fluorescent calcium indicators, 
the sole observation variable is the recorded intensity Ij (t) for cell j at a particular pixel reflecting 
some linear multiplier of the cytosolic calcium concentration at that point in the visual field, based 
on the dye loading in a cell, as well as the camera, microscope, and illumination setup. Defining 
— {^ji^)} with observation function F^, 

I,{t) = F,(c,(0,*,) = nc,it) + b + Ar{0,-f^) 

The parameter set #j = {71,6,7} represents the scaling, offset, and noise standard deviation of 
the observation function. The function A/'(0,7^) generates a normally distributed, random value 
with zero mean and 7 standard deviation. The noise term models the type of frame-to-frame 
variation typically seen in the amplification of the CCD signal prior to digitization. The size of 7 
is proportional to the magnitude of the noise, itself affected by camera type and gain settings. 

Since the framework is defined at cellular resolution, we are making the simplifying assumption 
that the recorded intensity represents the average intensity for the region of interest demarcating a 
specific cell j within a larger field. Additionally, if the camera records at a slower frequency than 
the transition dynamics f camera > where /camera is the camera recording frequency and At 

is the time incriment (c./. equation [5|, one is averaging intensity values for the duration that the 
camera shutter stays open. If this is the case then multiple sequential calcium concentration values 
would be averaged to produce a single intensity reading. 



3.4 GPU Implementation and Benchmarks 

The practical application of the theoretical framework both for simulation and, as will be described 
in a subsequent paper, for mapping the unknown functional connectivity of experimentally observed 
cellular networks necessitates its implementation in an appropriate computing environment. We 
have taken advantage of emerging high performance general purpose-graphics processing unit (GPU) 
parallel computing, although it can run as serial code on a normal central processing unit (CPU), 
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which we have also tested. Within the GPU environment, the code has been designed to run 
on nVIDIA graphics cards equipped with the CUDA interface (see http://www.nvidia.com/cuda). 
In this way, we can parallelize vertex dynamics, signaling dynaimcs, and observation integrations 
over many processor cores, achieving significant speedup over CPU or cluster computations. The 
framework and associated single cell dynamic and network connectivity models have been coded 
as compact MATLAB-callable libraries. All graphics user interface (GUI) and input /output (I/O) 
operations are handled through MATLAB and the code has been written in both MATLAB and 
plain C libraries that communicate through MATLAB. The libraries offer direct control over all 
parameters for all vertices. Using plain C language, any model that can be analytically described 
within the framework can be easily coded into a simulation library. From a practical experimental 
perspective GPU computation offers unprecedented scalability to larger systems with full access to 
all state variables and parameters, enabling rapid parallel simulation when the framework is run 
in the forward direction, and real-time dynamic mapping when the framework is applied to the 
inverse problem of mapping unknown functional connectivities of cellular neural networks. Speed 
and parallelization are critical for statistical simulation based identification methods such as particle 
filtering to operate in real time or near real time. 

Benchmarks in figure [4] show the relative speeds of the CPU and GPU implementations of 
two different dynamic cell models within the framework, an Izhikevitch model and a Hodgkin- 
Huxley model, simulated in 40 test networks (see section 3.5 below regarding the test networks). 
The parallel GPU implementation performed anywhere from 8 to 200 times faster than the same 
code executed on a single CPU thread. Performance was measured as a slowdown or speedup 
factor relative to real time. The Izhikevitch model had six state variables and a timestep of At 
1 millisecond. We simulated 10 seconds worth of data equivalent to 10,000 time steps. The more 
computationally intensive Hodgkin-Huxley model had eight state variables and a time step of 0.03 
milliseconds. One second worth of data was simulated for the HH model, corresponding to about 
33,333 time steps. Benchmarks were calculated as the dimensionless ratio of the actual period of 
time simulated (i.e. real time) to the amount of physical computational time it took the GPU 
or CPU to carry out the simulations. Because of the parallel implementation on the GPU, the 
performance falloff was much slower than on the CPU with increasing network size. In both 
CPU and GPU cases, performance decreased with increasing network density (total number of 
edges/vertices square), since in more dense networks more information is transferred between edges. 
This is more apparent in the random networks that have higher density than other network classes, 
producing relatively slower processing speeds. In general, the constant time step in the framework 
makes parallelization easy on GPU architectures, delivering network simulation performance near 
or even faster than real-time. For the HH model, the GPU was able to carry out the computations 
in essentially real time for all of the networks tested except for the largest random networks, which 
were about 10 tol5 times slower then real time for random networks with over a thousand vertices 
(upper left panel in figure |4|. For computationally simpler models such as the Izhikevitch model 
the GPU computations were actually faster then the period being simulated (i.e. faster than real 
time), ranging from roughly 10 to 20 times faster for most networks and about real time for the 
largest random networks (upper right panel in figure [4]). In contrast, CPU computations were 
always slower than real time, from 15 to 800 times slower for the HH model depending on the 
network class and size (lower left panel in figure |4]) and from just under near real time for small 10 
vertex networks to about 15 times slower for the largest random networks for the Izhikevitch model 
(lower right panel figure |4]). Additional GPU cards can further improve performance by splitting 
the task of advancing temporal cell dynamics; however transfer of information between cells is still 
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limited by memory and bus speeds, so dense networks will run slower than sparse networks. It 
is important to appreciate that the ability to carry out such forward simulations or to solve the 
inverse problem of mapping unknown functional connectivity topologies of networks in near real 
time or faster than real time using GPU computing is due to the mathematical construction of 
the framework and its efficient algorithmic implementation. It is impossible to do such network 
simulations or mappings in real time using biophysical compartmental simulation environments. 
The MATLAB and CUDA code for the framework are available for download from the authors' 
website (http://www.silva.ucsd.edu/Silva_Lab/Links.html). 



Hodgkin-Huxley Performance Relative to Real-Time Izhikevitch Simple Model Performance Relative to Real-Time 




10 30 100 300 1000 10 30 100 300 1000 

Netork Size (vertices) Netork Size (vertices) 



Figure 4: CPU and GPU benchmark results for framework simulations using Hodgkin-Huxley and 
Izhikevitch models of single cell dynamics. Two different network topologies for four different classes 
of networks were simulated (lattice, small world, scale free, and random), with network sizes (i.e. 
number of vertices) of 10, 30, 100, 300, and 1000 vertices simulated for each of the eight different 
networks. The top graphs show performance on a single nVIDIA Tesla C1060 GPU, expressed 
relative to the real time simulation period for each single cell dynamic model (see text), while the 
data for the bottom two graphs show the performance on a single core of the an Intel Core2Duo 
2.5GHz processor. Negative values represent a slowdown relative to real time, while positive values 
represent a speedup. All forty networks were tested on both. See section |3.5| below for details 
regarding the test networks and network classes. 



3.5 Standardized Tests for Connectivity Estimation 

Lastly, we propose a standardized basis test set to evaluate the effectiveness of mapping algorithms. 
A standard test set is well accepted in the field of non-linear optimization, providing a standard 
measure of different algorithms [53]. There are two benefits to having a standard set of networks 
to use for mapping. First, multiple algorithms can be evaluated against the same network and 
model, providing relative performance benchmarks. Second, data generated for a network using 
one dynamic model can be mapped using another dynamic model and comparisons can be made 
between the original topology and the mapped topology. This latter approach helps answer ques- 
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tions about model fitness, which are especially useful when trying to map data with multiple models 
or uncertainty in models. We emphasize that the test set we propose here should in no way be 
interpreted as implying that the full complexity and variability of real neural circuits and networks 
is captured or even described by the set. But we argue that any mathematical method or algorithm 
that claims to be able to deal with real cellular neural networks of any meaningful size (e.g. on 
the order of tens to hundreds of cells or larger) must at the very least be able to effectively and 
efficiently map functional networks derived by this test set, which offers a first order approximation 
to the dynamics and complexity displayed by such networks. 

Here we offer the foundation for such a test set: The location of vertices in physical space and 
their physical connections based on different connective classes. The choice of dynamic model, pa- 
rameters, and functional weights is left open and up to the discretion of the individual investigator, 
since they are specific to the network being studied and the mapping algorithms being designed, 
but can be directly implemented within the framework we have developed. Test networks vary in 
size from 10 to 1000 vertices, covering the range of cells that can be imaged simultaneously with 
fluorescence microscopy. At the small end of the scale, networks on the order of 10 vertices is about 
the limit of existing connectivity estimation methods [311 [H]- The upper end of the scale at a 1000 
vertices was chosen largely due to limits of computational power available at present. 

Each graph is composed of N interconnected vertices located in two or three dimensional physical 
space, with minimum distance constraints and other dynamic parameters as described below for 
each network class. The physical connectivity between vertices follow one of four different graph 
theoretical classifications, since there is no measurable and "real" network spatial geometry defined 
in this case; this is also discussed below. Vertices were positioned in geometric space randomly, 
but with a prescribed minimum distance between neighbors. We developed a simple algorithm to 
populate a physical space with TV cells or nodes: 

pick a random position for first cell XI 
set i=2 
while i<=N 

pick a random position Xi in space range 
set found=l 
for j=l to i-1 

if cell Xi is less than d units from cell Xj 
f ound=0 

break for loop 
end if 

end for 

if fomid==l 

accept Xi as next cell 
increment i 
end if 
end while 

The space range describes the dimensions of the space occupied by vertices in the network. For 
a two dimensional network this could be a square area of 500 x 500 distance units. Distance and 
position units are non-dimensional and can be scaled as needed. The minimum distance between 
vertices is a parameter, but can also be expressed as a function of the number of vertices N and 
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the dimension of the space the graph occupies: 

Here the minimum distance dmm is proportional to the length of the space L times the number 
of vertices N raised to the negative inverse of the dimension D. For example, a three dimensional 
space of length L = 100 could fit = 1000 vertices with minimum distance dmin — 100 x 1000^^/^. 
This implies a minimum distance of 10. While this is the absolute minimum distance with cubic 
packing, when placing vertices at random this minimum distance is reduced to allow for some 
variability in placement. Figure [5] shows examples of vertex placement in two and three dimensions 
for three different size networks. It is important to note that with this formula, the minimum 
distance can be prescribed for fractional dimensions, as may be the case for some neural tissues 
where cell arrangement is neither flat nor fully three dimensionally filling. Vertices are numbered 
from the center outwards, so mapping can be performed on a subset of vertices that interact with 
the complete network. This is a more difficult case but a more realistic scenario, as vertices would 
be receiving inputs from unobserved vertices in real cellular networks due to experimentally limited 
windows of observability. 

Once the vertices are placed in physical space, connections are made using established graph 
theoretic classes. We included lattice, small world, scale free, and random classes in the test set. 
For each connectivity class, we generated two networks of different edge densities, one with fewer 
and one with more edges. We chose these four classes because they represent the major graph the- 
oretic topologies, but of course any other algorithmically defined class can be used. The physical 
connectivity of a graph intuitively represents a constrained phase space on which dynamic signals 
propagate in both space and time (i.e. the functional connectivity topology) as determined by the 
network signaling framework and chosen model of single cell dynamics. Specifically, we considered 
the following classes and specific parameters for each: 

Lattice networks. This class of networks has only nearest neighbor connections with no long distance 
connections. The number of nearest local connections or total number of edges can be specified be- 
fore construction. In our case we limited connections for each vertex to its closest 3 and 8 neighbors. 

Small world networks. This is a modification of a lattice network which includes a specified prob- 
ability of long-range connections [3] . The probability of long range connections ranges between 
(lattice network) and 1 (random network), but typical values are around 10 percent, meaning that 
10 percent of all edges are randomly chosen. In our case we built networks of 5 and 15 percent 
probability of random re- wiring. 

Scale free networks. These networks follow a power law connection (edge) degree distribution, with 
many cells having few connections and few cells having many connections. When positional aspects 
are taken into account, scale free networks take on some small world properties and are essentially 
scale free geometric graphs called apoUonian networks [2]. 

Random networks. The study of random graphs extends all the way back to the original work of 
Erdos and Renyi. In a random graph, a specified number of edges are placed between randomly 
chosen vertices, without regards for vertex position, which has no meaning. We built random net- 
works of 10 and 20 percent densities, meaning about {N"^ * 10%) or (7V^ * 20%) number of edges 
where N is the number of vertices in the network. A 100% dense network connects every vertex to 
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Figure 5: Random placement of vertices with a minimum distance constraint. Examples of groups 
of 10, 100, and 1000 vertices placed in two (left column) and three (right column) dimensions 
are shown here. Vertices were colored randomly for clarity. The minimum distance is a function 
of both the number of vertices and the spatial dimension. For the same number of vertices, a 
two-dimensional network will have a smaller minimum distance and thus be packed tighter than a 
three-dimensional network. See text for details. 
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every other vertex. 

The different classes and densities are shown in figure [6) along with graph-theoretic statistics on 
connectivity and wiring lengths. The connections establish the physical connection between vertices 
or the edges along which functional connections are possible. The magnitude of the functional 
weights should be chosen based on the dynamic model and the specifics of the system studied. 
A mapping algorithm must identify the functional connections without knowledge of the physical 
connectivity class. Delays are defined according to the cartesian distances between connected 
vertices. Three-dimensional spaces generally have narrower and smaller distance distributions than 
two-dimensional packings. The formula presented in equation [25] is used to establish delays, with 
the speed parameter chosen arbitrarily based on the system being studied. The complete test set is 
80 networks, combinations of two dimensions, five sizes, and four connectivity classes, as outlined 
in table [l] Figure [7] shows an example of simulated calcium response raster plots for the networks 
shown in figure [6j A mapping algorithm should be able to identify the dynamics parameters and the 
functional connectivity of each test network, given a chosen single cell dynamic model. Ultimately, 
the only kind of measured experimental data available to any such algorithm would be imaged 
calcium responses such as those simulated in figure [7] or some equivalent data for another marker 
of functional cellular activity. These are the practical experimental constrains that any theoretical 
methods aimed at mapping functional activity in cellular neural circuits and networks with single 
cell resolution must face. The test networks are also available for download from the authors' 
website (http://www.silva.ucsd.edu/Silva_Lab/Links.html). 



Parameter 


Range 


Network Size N 


10, 30, 100, 300, 1000 


Geometric Dimension 


2D, 3D (may even be fractal, i.e. 2.5D) 


Connectivity Type 


lattice, small world, scale-free, random 


Edge Density 


low, high 



Table 1: Range of parameters specified by the test set. The current set is composed of eighty test 
networks span spanning ranges in network size, geometric dimension, connectivity type and edge 
density. The choice of dynamical cell model, parameters, functional weights, observation variable, 
user inputs and noise levels are left to the individual user. 



4 Discussion 

Within the study of networks, there are two opposing yet deeply interrelated processes: Simulation 
and estimation. Simulation of networks deals with the forward problem of making predictions 
using an established model and measured parameters and connectivities. The reverse problem, 
estimation or mapping, uses the result of actual collected data to infer, estimate, or map the 
parameters of a model or, for networks, functional connectivity. The framework we introduce 
here can be used both for simulating signal propagation in physically realistic networks, and for 
the reverse process of estimating or mapping unknown functional connectivities of networks. The 
framework is bounded by a set of rules and constraints imposed by the experimental reality of 
cellular neurobiological methods: Complex non-linear dynamics, limited observability, noise and 
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Figure 6: Test network topologies for a 100-vertex two dimensional network. For each of the 
four topologies included in the test set two edge densities were considered, one with fewer edges 
or randomness (green), and a more complex one with more edges or randomness (blue; left two 
columns). Parameter values for each network topology are shown below their respective graphs. The 
two right-most columns plot the degree distribution (third column from the left) and the geometric 
length distribution (fourth column from the left) for each of the edge densities for each class. The 
total number of edges for each network is the area under the edge distribution plots, while the total 
physical wiring length of the network is the area under the geometric length distribution plots. 
The higher edge density networks (colored in blue) have both higher total number of edges and 
total wiring length than the lower edge density networks (colored in green). For a constant signal 
propagation speed the wiring length distribution represents the delay time distribution between 
vertices. 
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Time (s) Time (s) 

Figure 7: Computed (i.e. simulated) dynamics of experimentally observable calcium signaling for 
each of the eight networks shown in figure [6j arranged in the same order. Each network was 
simulated for four seconds using the Izhikevitch simple model, with a constant current input to 
generate activity. The gray scale represents calcium activity with lighter shades representing the 
approximate instantaneous spike rate of individual neurons in the network. 
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uncertainty, and experimental control. It is designed around current observation and experimental 
capabilities, which are shifting from single neuron multiple trials, to multiple neuron single trial 
experiments |36j . 

Within the framework, the mathematical construction for the dynamic model describes the time 
course of each vertex. A general state transition representation encompasses different model types, 
from ordinary differential equations to state machines and Markov models, to simplified neuronal 
models like LIF and Izhikevitch. The choice of model will certainly affect the estimation of the 
dynamic parameters from the collected data; how much the single-cell dynamic model affects the 
estimation of functional weights is is still an open question. The proposed test set should help 
address this by simulating artificial data with one model and estimating with another. 

Mapping a complete functional topology is ultimately a reverse process, and will involve some 
combination of estimation, filtering, and optimization. While some approaches exist for estimating 
parameters and dynamics of single neurons [13] or small groups of neurons [311 [171 [H] j mapping 
large networks within biologically realistic constraints remains a challenge, and we are still a long 
way from establishing a complete functional connectivity map of even simple processes and tasks. 
Indeed, from a neurophysiological perspective, it is not even entirely clear 'what' we should be 
mapping or how to properly interpret such data from the perspective of deciphering the neural 
code. The proposed framework attempts to unify both theoretical and practical considerations as 
an "open standard" for the development of large scale functional topology reconstruction algorithms. 

Since the goal of mapping is to identify both the dynamic parameters of individual vertices as 
well as the connectivity between vertices, a well-designed input control should be used to make the 
observation as informative as possible, provided it does not alter the parameters and connectivities 
of interest. Experimentally, input control can take on many forms. The dynamics of individual 
cells can be perturbed using methods such as optogenetics, pharmacologically using appropriate 
agonists and antagonists, and electrophysiology. For single cells there are different input functions 
that can be used and there are a few approaches describing input function design to extract the 
most amount of information [521 H] • At the network level, the set of input functions for each cell 
must be designed in parallel and coordinated with observed activity in order to provide the most 
amount of information to the mapping algorithm. 
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